clear
rng default

all_years=1940:1:1967;
n_years=size(all_years,2);
ntypes=4;

C_21_m_Gumbel=zeros(n_years,1);
C_32_m_Gumbel=zeros(n_years,1);
C_43_m_Gumbel=zeros(n_years,1);
C_21_m_Gumbel_original=zeros(n_years,1);
C_32_m_Gumbel_original=zeros(n_years,1);
C_43_m_Gumbel_original=zeros(n_years,1);


load Cm_R1
C_21_m_R1=C_21_m(1:n_years,:);
C_32_m_R1=C_32_m(1:n_years,:);
C_43_m_R1=C_43_m(1:n_years,:);

load Cm_R2
C_21_m_R2=C_21_m(1:n_years,:);
C_32_m_R2=C_32_m(1:n_years,:);
C_43_m_R2=C_43_m(1:n_years,:);




for j=1:n_years
    
year=all_years(j);

% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond=[condprobsmen(:,2:end) condprobsmen(:,1)]; %4x5
%Note: PYX_cond(x,y)=Probability to marry woman y born in any year conditional on being man x born in 1940

% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year+1) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year+1) '_140']);
PXY_cond=[condprobswomen(:,2:end) condprobswomen(:,1)];  %4x5
%Note: PXY_cond(y,x)=Probability to marry man x born in any year conditional on being woman y born in 1941

% Probability distribution of X
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX=probsmen_temp.(['probsmen' num2str(year) '_140']);
%PX(x)=Probability of being man x conditional on being born in 1940

% Probability distribution of Y
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year+1) '_140']);
PY=probswomen_temp.(['probswomen' num2str(year+1) '_140']);
%PY(y)=Probability of being woman y conditional on being born in 1941

[U_Gumbel, ~, ~, ~,~,~]=param_Gumbel_closedform(PYX_cond, PXY_cond, ntypes, ntypes);

Mmarital_systematic = zeros(1,ntypes-1);
for h=1:ntypes-1
    Mmarital_systematic(j,h)=sum(PYX_cond(h+1,1:end-1).*U_Gumbel(h+1,:),2)-sum(PYX_cond(h,1:end-1).*U_Gumbel(h,:),2);
end

Mmarital = zeros(1,ntypes-1);
for h=1:ntypes-1
    Mmarital(j,h)=(log(exp(U_Gumbel(h+1,1))+exp(U_Gumbel(h+1,2))+exp(U_Gumbel(h+1,3))+exp(U_Gumbel(h+1,4))+1)+ 0.5772)- ...
                  (log(exp(U_Gumbel(h,1))+exp(U_Gumbel(h,2))+exp(U_Gumbel(h,3))+exp(U_Gumbel(h,4))+1)+ 0.5772);
end


C_21_m_Gumbel(j)=Mmarital_systematic(j,1);
C_32_m_Gumbel(j)=Mmarital_systematic(j,2);
C_43_m_Gumbel(j)=Mmarital_systematic(j,3);

C_21_m_Gumbel_original(j)=Mmarital(j,1);
C_32_m_Gumbel_original(j)=Mmarital(j,2);
C_43_m_Gumbel_original(j)=Mmarital(j,3);

end


%% Figure 5 in the paper
figure
plot(all_years, C_21_m_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_21_m_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{21}(U) and \\Delta_{21}(U)'),'FontSize',15)
saveas(gcf, 'C_21_m_original.jpg')


figure
plot(all_years, C_32_m_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_32_m_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{32}(U) and \\Delta_{32}(U)'),'FontSize',15)
saveas(gcf, 'C_32_m_original.jpg')


figure
plot(all_years, C_43_m_Gumbel, 'b', 'LineWidth',3)
hold on
plot(all_years, C_43_m_Gumbel_original, 'k', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
box on
xlim([1940 1967])
ylim([ -1 1 ])
set(gca,'ytick', -1:0.5:1 , 'FontSize',15)
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{43}(U) and \\Delta_{43}(U)'),'FontSize',15)
saveas(gcf, 'C_43_m_original.jpg')



%% Figure 6 in the paper

C_21_m_min=C_21_m_R1(:,1);
C_21_m_min(C_21_m_min==-Inf)=[];
if isempty(C_21_m_min)
    C_21_m_min=min(C_21_m_Gumbel)-2;
else
    C_21_m_min=min(C_21_m_min)-2;
end
C_21_m_max=C_21_m_R1(:,2);
C_21_m_max(C_21_m_max==Inf)=[];
if isempty(C_21_m_max)
    C_21_m_max=max(C_21_m_Gumbel)+2;
else
    C_21_m_max=max(C_21_m_max)+2;
end


%%
C_32_m_min=C_32_m_R1(:,1);
C_32_m_min(C_32_m_min==-Inf)=[];
if isempty(C_32_m_min)
    C_32_m_min=min(C_32_m_Gumbel)-2;
else
    C_32_m_min=min(C_32_m_min)-2;
end
C_32_m_max=C_32_m_R1(:,2);
C_32_m_max(C_32_m_max==Inf)=[];
if isempty(C_32_m_max)
    C_32_m_max=max(C_32_m_Gumbel)+2;
else
    C_32_m_max=max(C_32_m_max)+2;
end


%%
C_43_m_min=C_43_m_R1(:,1);
C_43_m_min(C_43_m_min==-Inf)=[];
if isempty(C_43_m_min)
    C_43_m_min=min(C_43_m_Gumbel)-2;
else
    C_43_m_min=min(C_43_m_min)-2;
end
C_43_m_max=C_43_m_R1(:,2);
C_43_m_max(C_43_m_max==Inf)=[];
if isempty(C_43_m_max)
    C_43_m_max=max(C_43_m_Gumbel)+2;
else
    C_43_m_max=max(C_43_m_max)+2;
end

y_axis_min=round(min([C_21_m_min C_32_m_min C_43_m_min]));
y_axis_max=round(max([C_21_m_max C_32_m_max C_43_m_max]));



%% 21

for h=1:n_years
    if  C_21_m_R1(h,1)==-Inf
        C_21_m_R1(h,1)=y_axis_min;
    end
    if  C_21_m_R2(h,1)==-Inf
        C_21_m_R2(h,1)=y_axis_min;
    end
end


for h=1:n_years
    if  C_21_m_R1(h,2)==Inf
        C_21_m_R1(h,2)=y_axis_max;
    end
    if  C_21_m_R2(h,2)==Inf
        C_21_m_R2(h,2)=y_axis_max;
    end
end




figure
plot(all_years, C_21_m_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=C_21_m_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_21_m_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_21_m_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_21_m_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_21_m_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
%ytickformat('%,.1f')
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{21}(U)'),'FontSize',15)
saveas(gcf, 'C_21_m.jpg')


%% 32

for h=1:n_years
    if  C_32_m_R1(h,2)==Inf
        C_32_m_R1(h,2)=y_axis_max;
    end
    if  C_32_m_R2(h,2)==Inf
        C_32_m_R2(h,2)=y_axis_max;
    end
end

for h=1:n_years
    if  C_32_m_R1(h,1)==-Inf
        C_32_m_R1(h,1)=y_axis_min;
    end
    if  C_32_m_R2(h,1)==-Inf
        C_32_m_R2(h,1)=y_axis_min;
    end
end


figure
plot(all_years, C_32_m_R1(:,2), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=C_32_m_R1(:,1).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_32_m_R1(:,2).', fliplr(upline)];
fill(polyX, polyY, [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_32_m_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_32_m_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_32_m_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
%ytickformat('%,.1f')
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{32}(U)'),'FontSize',15)
saveas(gcf, 'C_32_m.jpg')


%% 43

for h=1:n_years
    if  C_43_m_R1(h,1)==-Inf
        C_43_m_R1(h,1)=y_axis_min;
    end
    if  C_43_m_R2(h,1)==-Inf
        C_43_m_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  C_43_m_R1(h,2)==Inf
        C_43_m_R1(h,2)=y_axis_max;
    end
    if  C_43_m_R2(h,2)==Inf
        C_43_m_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, C_43_m_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=C_43_m_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_43_m_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=C_43_m_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [C_43_m_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, C_43_m_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([y_axis_min y_axis_max])
set(gca,'ytick',y_axis_min:3:y_axis_max, 'FontSize',15)
%ytickformat('%,.1f')
xlabel('Cohort','FontSize',15) 
ylabel(sprintf('C_{43}(U)'),'FontSize',15)
saveas(gcf, 'C_43_m.jpg')


close all